LoadLibrary <- function(){
  library(plyr)
  library(zoo)
  library(xts)
  library(dygraphs)
  library(lubridate)
  library(colorspace)
  library(ggplot2)
  library(RColorBrewer)
  library(forecast)
  library(data.table)
  library(MASS)
  library(gcookbook)
  library(scales)
  library(gdata)
  library(plm) # Panel Data Analysis
  library(stargazer)
}


DefineTableNA <- function(x){
  x[x == -9999] <- NA 
}


GetYearMonth <- function(x) {
  x$year <- year(x$date)
  x$month <- months(x$date, abbreviate = T)
  x$month <- match(x$month, month.abb)
  x <- x[order(x$year, x$month),] 
  x
}


OrderByStateYear <- function(x){
  x <- x[order(x$state, x$year),]
}


SelectYear1995_2014 <- function(d){
  d <- subset(d, year >= 1995 & year<= 2014)
}


ConvertToCelcius <- function(x){
  x <- x/10
}

ConvertToMM <- function(x){
  x <- x/10
}

RoundTo_num_DecimalPoints <- function(x, num){
  x <- round(x,num)
}


PlotFittedYield <- function(reg, title){
  y_fitted <- reg$model[ ,1] - reg$residuals
  y_fitted <- as.numeric(y_fitted)
  d <- data.frame(reg$model[ ,c("grain_yi","time_trend")], y_fitted)
  names(d) <- c("y", "year", "y_fitted")
  dforplot <- melt(d, id.vars="year")
  dforplot[,1] <- as.numeric(dforplot[,1])
  dforplot[,3] <- as.numeric(dforplot[,3])
  
  ggplot(dforplot, aes(x=year, y=value, col=variable)) + 
    geom_line(size=0.5) +
    geom_point(size=0.5) +
    geom_smooth(size=0.5) + # Loess smooth; se = FALSE can be added to remove bandwidth shadow
    scale_x_continuous(breaks=seq(1995,2014,1)) +
    theme(axis.text.x = element_text(angle = 30)) +
    theme(axis.line = element_line(colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank() ) +
    ggtitle(paste0(title,", 1995-2014")) + 
    xlab("Year") + 
    ylab("Yield (bu/acre)")
}


PlotFittedProduction <- function(reg, title){
  y_fitted <- reg$model[ ,1] - reg$residuals
  y_fitted <- as.numeric(y_fitted)
  d <- data.frame(reg$model[ ,c("grain_prod","time_trend")], y_fitted)
  names(d) <- c("y", "year", "y_fitted")
  dforplot <- melt(d, id.vars="year")
  dforplot[,1] <- as.numeric(dforplot[,1])
  dforplot[,3] <- as.numeric(dforplot[,3])
  
  ggplot(dforplot, aes(x=year, y=value, col=variable)) + 
    geom_line(size=0.5) +
    geom_point(size=0.5) +
    geom_smooth(size=0.5) + # Loess smooth; se = FALSE can be added to remove bandwidth shadow
    scale_x_continuous(breaks=seq(1995,2014,1)) +
    theme(axis.text.x = element_text(angle = 30)) +
    theme(axis.line = element_line(colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank() ) +
    ggtitle(paste0(title,", 1995-2014")) + 
    xlab("Year") + 
    ylab("Production (bu)")
}
